Last updated: 2025-11-16

Checks: 5 2

Knit directory: PIPAC_spatial/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: uncommitted changes

The R Markdown is untracked by Git. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20240917)

The command set.seed(20240917) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: absolute

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute relative
/home/hnatri/PIPAC_spatial/ .
/home/hnatri/PIPAC_spatial/code/PIPAC_colors_themes.R code/PIPAC_colors_themes.R
/home/hnatri/PIPAC_spatial/code/plot_functions.R code/plot_functions.R
/home/hnatri/PIPAC_spatial/cell_nonimmune_cluster_marker_annotations.tsv cell_nonimmune_cluster_marker_annotations.tsv

Repository version: 2221fdb

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 2221fdb. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 


Ignored files:
    Ignored:    analysis/figure/
    Ignored:    cell_immune_cluster_marker_annotations.tsv
    Ignored:    cell_nonimmune_cluster_marker_annotations.tsv
    Ignored:    celltype_markers.tsv
    Ignored:    code/.Rhistory
    Ignored:    code/Proximity_analysis/.DS_Store
    Ignored:    code/Proximity_analysis/data/
    Ignored:    code/Proximity_analysis/output/
    Ignored:    code/RSC_latest_EDM_2025-08-06/
    Ignored:    code/celltype_proximity_tumor.Rout
    Ignored:    immune_cluster_marker_annotations_2ndpass.tsv
    Ignored:    immune_prelim_annot_topmarkers.tsv
    Ignored:    nonimmune_cluster_marker_annotations_2ndpass.tsv
    Ignored:    nonimmune_prelim_annot_topmarkers.tsv

Untracked files:
    Untracked:  analysis/finalize_annotations.Rmd
    Untracked:  analysis/immune_annotation_cell.Rmd
    Untracked:  analysis/nonimmune_annotation_cell.Rmd
    Untracked:  code/.RDataTmp
    Untracked:  code/.RDataTmp1
    Untracked:  code/celltype_proximity_09232025.R
    Untracked:  code/celltype_proximity_tumor.R
    Untracked:  code/denoist_chunk.R
    Untracked:  code/plot_informative_cluster_markers.R
    Untracked:  code/slurm.26717417.err
    Untracked:  code/slurm.26717417.out
    Untracked:  enrichment_posttreatment_scores_compiled.tsv
    Untracked:  enrichment_pretreatment_scores_compiled.tsv
    Untracked:  enrichment_tumor_scores_compiled.tsv
    Untracked:  output/enrichment_scores_compiled.tsv
    Untracked:  output/proximity_compiled.tsv
    Untracked:  tumor_posttreatment_proximity_compiled.tsv
    Untracked:  tumor_pretreatment_proximity_compiled.tsv
    Untracked:  tumor_proximity_compiled.tsv

Unstaged changes:
    Deleted:    Rplots.pdf
    Modified:   analysis/Xenium_preprocess_wholecell.Rmd
    Modified:   analysis/add_metadata.Rmd
    Modified:   analysis/annotation_cell.Rmd
    Modified:   analysis/arm3_prepost_plots_DEGs.Rmd
    Modified:   analysis/compare_cellular_nuclear.Rmd
    Modified:   analysis/compare_cellular_nuclear_denoist_annot.Rmd
    Modified:   analysis/index.Rmd
    Modified:   analysis/pca_variance_decomp.Rmd
    Modified:   analysis/proximity.Rmd
    Modified:   analysis/splitting_samples.Rmd
    Modified:   code/PIPAC_colors_themes.R
    Modified:   code/anndata_to_seurat.R
    Modified:   code/denoist.R
    Modified:   code/pairwise_proximity.R
    Modified:   code/pairwise_proximity_chunk.R
    Modified:   code/parse_denoist_res.R
    Modified:   code/plot_functions.R
    Modified:   code/run_rscript.sh
    Modified:   code/seurat_to_anndata.R
    Modified:   code/update_metadata.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

There are no past versions. Publish this analysis with wflow_publish() to start tracking its development.

Packages and environment variables

suppressPackageStartupMessages({
  library(workflowr)
  library(arrow)
  library(dplyr)
  library(Seurat)
  library(SeuratObject)
  library(SeuratDisk)
  library(tidyverse)
  library(tibble)
  library(ggplot2)
  library(ggpubr)
  library(ggrepel)
  library(googlesheets4)
  library(workflowr)})

Environment variables and helper functions

setwd("/home/hnatri/PIPAC_spatial/")
set.seed(9999)
options(scipen = 99999)
options(ggrepel.max.overlaps = Inf)

source("/home/hnatri/PIPAC_spatial/code/PIPAC_colors_themes.R")
source("/home/hnatri/PIPAC_spatial/code/plot_functions.R")

recluster_seurat <- function(seurat_subset, umap_pcs = 10){
    seurat_subset <- SCTransform(seurat_subset,
                                 vars.to.regress = c("nCount_RNA", "nFeature_RNA"),
                                 verbose = F)
      
    seurat_subset <- RunPCA(seurat_subset,
                            npcs = 50,
                            assay = "SCT",
                            #features = features,
                            approx = F,
                            verbose = F)
    
    integrated_data <- RunUMAP(seurat_subset,
                               reduction = "pca",
                               reduction.name = "umap",
                               dims = 1:umap_pcs,
                               return.model = TRUE,
                            verbose = F)

    integrated_data <- FindNeighbors(integrated_data,
                                     reduction = "pca",
                                     dims = 1:umap_pcs,
                                     graph.name = c("SCTnn",
                                                    "SCTsnn"),
                            verbose = F)

    integrated_data <- FindClusters(integrated_data,
                                    resolution = c(0.1,0.2,0.3,0.5,0.8,1),
                                    graph.name = "SCTsnn",
                            verbose = F)
    
    return(integrated_data)
}

Marker information and metadata

# Cell type markers by group
gs4_deauth()
markers_by_group  <- gs4_get("https://docs.google.com/spreadsheets/d/1w0wrL-5KwNEWi_VpmriN7YmfwpBPlEYX0MKlwI2ZQu8/edit?usp=sharing")
markers_by_group <- read_sheet(markers_by_group, sheet = "Markers by group")

metadata  <- gs4_get("https://docs.google.com/spreadsheets/d/1sXXwOreLxjMSUoPt79c6jmaQpluWkaxA5P5HfDsed3I/edit?usp=sharing")
markers <- read_sheet(metadata, sheet = "Markers")
first_pass_annot <- read_sheet(metadata, sheet = "Cell based annotation, non-immune")

Import data

seurat_data <- readRDS("/tgen_labs/banovich/PIPAC/Seurat/cell_nonimmune_clustered_NN30_PC20_Seurat.rds")

DefaultAssay(seurat_data) <- "denoist_RNA"
seurat_data <- NormalizeData(seurat_data)

DimPlot(seurat_data,
        group.by = "leiden_1.5",
        cols = pipac_cell_nonimmune_cluster_col,
        reduction = "umap",
        raster = T,
        label = T) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend()

QC metrics for each cluster

VlnPlot(seurat_data,
        features = c("percent.mt", "nCount_cell_RNA", "nFeature_cell_RNA"),
        group.by = "leiden_1.5",
        ncol = 1,
        pt.size = 0) &
  theme_bw() &
  NoLegend()

seurat_data$leiden_1.5chr <- as.character(seurat_data$leiden_1.5)
seurat_data$leiden_1.5chr <- factor(seurat_data$leiden_1.5chr,
                                    levels = as.character(sort(unique(seurat_data$leiden_1.5))))
DimPlot(seurat_data,
        group.by = "leiden_1.5chr",
        split.by = "leiden_1.5chr",
        cols = pipac_cell_nonimmune_cluster_col,
        reduction = "umap",
        raster = T,
        #label = T,
        ncol = 6) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend()

Top cluster markers

Idents(seurat_data) <- seurat_data$leiden_1.5
cluster_markers <- FindAllMarkers(seurat_data,
                                  return.thresh = 0.01,
                                  logfc.threshold = 0.5,
                                  min.pct = 0.20,
                                  verbose = T)

table(cluster_markers$cluster)

  7   2  26  15  16   3  21  29  31  12  35  32  25   8  22   6  10  36  24  13 
 37  43  23  75  25  32  30  50  21  21 140 122  82  75  28  35  25  54  46  38 
  4  30  33   0  20  23  18  19   9  28  34  17   1  11  27  14   5 
 44  25 142  34  26  42  86  55  54  27  23  37  33  38  47  52  27 
#hist(cluster_markers$avg_log2FC, main = "", xlab = "avg_log2FC", breaks = 100)
#hist(cluster_markers$p_val, main = "", xlab = "p_val", breaks = 100)
#hist(cluster_markers$p_val_adj, main = "", xlab = "p_val_adj", breaks = 100)

top_cluster_markers <- cluster_markers %>%
  arrange(dplyr::desc(avg_log2FC)) %>%
  group_by(cluster) %>%
  dplyr::slice(1:5)
plot_features <- c("PTPRC",
                   "CD3D", "CD3E", "CD4", "CD8A", # T cells
                   "STAT4", "STAT3", "TIGIT", "GZMB",
                   "SELL", "CD19", # B cells
                   "CD68", "CD44", "MARCO", "APOE", # Macrophages
                   "C1QB", "C1QBP",
                   "MUC5AC", "NOTCH3", "RGS5", "MS4A1", "PGA5", "PLVAP", # Lineage markers
                   "FN1", "DCN", "LUM", # Fibroblasts
                   "EGR3", "TP53", "JUN", "KIT", # Tumor
                   "SOX9", "RNF43", "EPCAM")

seurat_data$leiden_1.5 <- factor(seurat_data$leiden_1.5,
                                 levels = sort(unique(seurat_data$leiden_1.5)))
#Idents(seurat_data) <- seurat_data$leiden_1.5

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = plot_features,
        cols = c("gray98", "tomato3")) +
  RotatedAxis()

FeaturePlot(seurat_data,
            features = plot_features,
            cols = c("gray98", "tomato3"),
            reduction = "umap",
            raster = T,
            ncol = 5) &
  coord_fixed(ratio = 1) &
  theme_bw() &
  NoLegend()

create_dotplot_heatmap(seurat_object = seurat_data,
                       plot_features = unique(top_cluster_markers$gene),
                       group_var = "leiden_1.5",
                       group_colors = pipac_cell_nonimmune_cluster_col,
                       column_title = "",
                       row_km = 5,
                       col_km = 5,
                       row.order = NULL,
                       col.order = NULL)

Comparing to nuclear annotations

create_barplot(seurat_data,
               plot_var = "Annotation",
               group_var = "leiden_1.5",
               plot_levels = sort(unique(seurat_data$Annotation)),
               group_levels = sort(unique(seurat_data$leiden_1.5)),
               plot_colors = pipac_celltype_col,
               var_names = c("Annotation", "leiden_1.5"),
               legend_title = "")

create_barplot(seurat_data,
               group_var = "Annotation",
               plot_var = "leiden_1.5",
               group_levels = sort(unique(seurat_data$Annotation)),
               plot_levels = sort(unique(seurat_data$leiden_1.5)),
               plot_colors = pipac_cell_nonimmune_cluster_col,
               var_names = c("leiden_1.5", "Annotation"),
               legend_title = "")

Saving top markers and annotations

output_cluster_markers <- cluster_markers %>%
  arrange(dplyr::desc(avg_log2FC)) %>%
  group_by(cluster) %>%
  dplyr::slice(1:30)

output_cluster_markers <- merge(output_cluster_markers, markers, by.x = "gene", by.y = "Gene")

write.table(output_cluster_markers, "/home/hnatri/PIPAC_spatial/cell_nonimmune_cluster_marker_annotations.tsv",
            quote = F, row.names = F, sep = "\t")

Plotting cell type markers by group

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Cell-cell crosstalk"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Cell-cell crosstalk")

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Macrophage signaling"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Macrophage signaling")

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "DNA damage and cell faith"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("DNA damage and cell faith")

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Functional markers"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Functional markers")

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Mesothelial"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Mesothelial")

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Macrophages"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Macrophages")
DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "CAFs"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("CAFs")

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = grep("COL", rownames(seurat_data), value = T),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Collagen")

Cell cycle markers

s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes

DotPlot(seurat_data,
        group.by = "leiden_1.5",
        features = c(s.genes, g2m.genes),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("S/G2M genes")

Further splitting and clustering non-malignant cells

tumor_clusters <- c(7, 8, 25, 32, 33, 35)
stromal_seurat_data <- subset(seurat_data, subset = leiden_1.5 %in% tumor_clusters, invert = T)

#stromal_seurat_data <- recluster_seurat(stromal_seurat_data, vars_to_regress = NULL)

#saveRDS(stromal_seurat_data, "/scratch/hnatri/PIPAC/stromal_seurat_data.rds")
stromal_seurat_data <- readRDS("/tgen_labs/banovich/PIPAC/Seurat/stromal_clustered_NN20_PC20_Seurat.rds")
DefaultAssay(stromal_seurat_data) <- "denoist_RNA"
stromal_seurat_data <- NormalizeData(stromal_seurat_data)

DimPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        cols = pipac_cell_nonimmune_cluster_col,
        reduction = "umap",
        raster = T,
        label = T) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend()

Marker expression in the non-malignant subset

Idents(stromal_seurat_data) <- stromal_seurat_data$leiden_1.0
cluster_markers <- FindAllMarkers(stromal_seurat_data,
                                  return.thresh = 0.01,
                                  logfc.threshold = 0.5,
                                  min.pct = 0.20,
                                  verbose = T)

top_cluster_markers <- cluster_markers %>%
  arrange(dplyr::desc(avg_log2FC)) %>%
  group_by(cluster) %>%
  dplyr::slice(1:5)

Plotting lineage markers

plot_features <- c("PTPRC",
                   "CD3D", "CD3E", "CD4", "CD8A", # T cells
                   "STAT4", "STAT3", "TIGIT", "GZMB",
                   "SELL", "CD19", # B cells
                   "CD68", "CD44", "MARCO", "APOE", # Macrophages
                   "C1QB", "C1QBP", "CD163",
                   "SPP1", "VCAN",
                   "MUC5AC", "NOTCH3", "RGS5", "MS4A1", # Lineage markers
                   "PGA5", "PLVAP", "CCL21", "LYVE",
                   "SPARCL1", # Endo/peri
                   "FN1", "DCN", "LUM", # Fibroblasts
                   "COL1A2", "PDPN", "ACTA2", "LMNA",
                   "EGR3", "TP53", "JUN", "KIT", # Tumor
                   "SOX9", "RNF43", "EPCAM", "PECAM", "CEACAM")

stromal_seurat_data$leiden_1.0 <- factor(stromal_seurat_data$leiden_1.0,
                                 levels = sort(unique(stromal_seurat_data$leiden_1.0)))
Idents(stromal_seurat_data) <- stromal_seurat_data$leiden_1.0

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = plot_features,
        cols = c("gray98", "tomato3")) +
  RotatedAxis()

FeaturePlot(stromal_seurat_data,
            features = plot_features,
            cols = c("gray98", "tomato3"),
            reduction = "umap",
            raster = T,
            ncol = 5) &
  coord_fixed(ratio = 1) &
  theme_bw() &
  NoLegend()

Plotting top markers

create_dotplot_heatmap(seurat_object = stromal_seurat_data,
                       plot_features = unique(top_cluster_markers$gene),
                       group_var = "leiden_1.0",
                       group_colors = pipac_cell_nonimmune_cluster_col,
                       column_title = "",
                       row_km = 5,
                       col_km = 5,
                       row.order = NULL,
                       col.order = NULL)

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Cell-cell crosstalk"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Cell-cell crosstalk")

DotPlot(seurat_data,
        group.by = "leiden_1.0",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "DNA damage and cell faith"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("DNA damage and cell faith")

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Functional markers"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Functional markers")

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "Mesothelial"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Mesothelial")

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = unique(markers_by_group[which(markers_by_group$annotation1 == "CAFs"),]$feature),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("CAFs")

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = grep("COL", rownames(seurat_data), value = T),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("Collagen")

Cell cycle markers

s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes

DotPlot(stromal_seurat_data,
        group.by = "leiden_1.0",
        features = c(s.genes, g2m.genes),
        cols = c("azure", "tomato3")) +
  RotatedAxis() +
  ggtitle("S/G2M genes")

Cluster overlap with previous annotations

pipac_celltype_col <- c(pipac_celltype_col, c("NA" = "lightgray"))
stromal_seurat_data$Annotation <- ifelse(is.na(stromal_seurat_data$Annotation), "NA", stromal_seurat_data$Annotation)

create_barplot(seurat_data,
               group_var = "Annotation",
               plot_var = "leiden_0.5",
               group_levels = sort(unique(seurat_data$Annotation)),
               plot_levels = sort(unique(seurat_data$leiden_1.0)),
               plot_colors = pipac_cluster_20_col,
               var_names = c("leiden_0.5", "Annotation"),
               legend_title = "")

create_barplot(seurat_data,
               plot_var = "Annotation",
               group_var = "leiden_1.0",
               plot_levels = sort(unique(seurat_data$Annotation)),
               group_levels = sort(unique(seurat_data$leiden_1.0)),
               plot_colors = pipac_celltype_col,
               var_names = c("Annotation", "leiden_0.5"),
               legend_title = "")

create_barplot(stromal_seurat_data,
               group_var = "Annotation",
               plot_var = "leiden_0.5",
               group_levels = sort(unique(stromal_seurat_data$Annotation)),
               plot_levels = sort(unique(stromal_seurat_data$leiden_1.0)),
               plot_colors = pipac_cluster_20_col,
               var_names = c("leiden_0.5", "Annotation"),
               legend_title = "")

create_barplot(stromal_seurat_data,
               plot_var = "Annotation",
               group_var = "leiden_1.0",
               plot_levels = sort(unique(stromal_seurat_data$Annotation)),
               group_levels = sort(unique(stromal_seurat_data$leiden_1.0)),
               plot_colors = pipac_celltype_col,
               var_names = c("Annotation", "leiden_0.5"),
               legend_title = "")

Reclustering difficult clusters from the non-malignant subset

stromal_recluster_subset <- subset(stromal_seurat_data, subset = leiden_1.0 %in% c(0, 15, 23, 26))
stromal_recluster_subset <- recluster_seurat(stromal_recluster_subset, umap_pcs = 10)

DimPlot(stromal_recluster_subset,
        group.by = "SCTsnn_res.0.8",
        label = T) +
  NoLegend()

stromal_recluster_subset$SCTsnn_res.0.8 <- factor(stromal_recluster_subset$SCTsnn_res.0.8,
                                 levels = sort(unique(stromal_recluster_subset$SCTsnn_res.0.8)))
Idents(stromal_recluster_subset) <- stromal_recluster_subset$SCTsnn_res.0.8

DotPlot(stromal_recluster_subset,
        #group.by = "leiden_1.0",
        features = plot_features,
        cols = c("gray89", "tomato3")) +
  RotatedAxis()

Idents(stromal_recluster_subset) <- stromal_recluster_subset$SCTsnn_res.0.8
stromal_cluster_markers <- FindAllMarkers(stromal_recluster_subset,
                                          return.thresh = 0.01,
                                          logfc.threshold = 0.5,
                                          min.pct = 0.20,
                                          verbose = T)

table(stromal_cluster_markers$cluster)

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
128  73  37  69  97  34  75  66 121  86  39  73 119 101  73 104  93  71  60  79 
 20 
113 
stromal_top_cluster_markers <- stromal_cluster_markers %>%
  arrange(dplyr::desc(avg_log2FC)) %>%
  group_by(cluster) %>%
  dplyr::slice(1:5)

create_dotplot_heatmap(seurat_object = stromal_recluster_subset,
                       plot_features = unique(stromal_top_cluster_markers$gene),
                       group_var = "SCTsnn_res.0.8",
                       group_colors = pipac_cell_nonimmune_cluster_col,
                       column_title = "",
                       row_km = 5,
                       col_km = 5,
                       row.order = NULL,
                       col.order = NULL)

Different or same cells expressing multiple lineage markers?

Cluster 10

FeaturePlot(subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 10),
            features = c("DCN", "PLVAP"),
            blend = T) &
  theme_bw()

cl10 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 10)
cl10 <- recluster_seurat(cl10, umap_pcs = 10)

FeaturePlot(cl10,
            features = c("DCN", "PLVAP"),
            blend = T) &
  theme_bw()

FeaturePlot(cl10,
            features = c("DCN", "PLVAP")) &
  theme_bw()

ncol(cl10)
[1] 1102
table(cl10$SCTsnn_res.0.8)

  0   1   2   3   4   5   6 
281 183 183 146 135 106  68 
DimPlot(cl10,
        group.by = "SCTsnn_res.0.8",
        reduction = "umap",
        label = T) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend() +
  ggtitle("Cluster 10")

cl10$SCTsnn_res.0.8 <- factor(cl10$SCTsnn_res.0.8,
                              levels = sort(unique(cl10$SCTsnn_res.0.8)))
Idents(cl10) <- cl10$SCTsnn_res.0.8

DotPlot(cl10,
        group.by = "SCTsnn_res.0.8",
        features = plot_features,
        cols = c("gray98", "tomato3")) +
  RotatedAxis()

Idents(cl10) <- cl10$SCTsnn_res.0.8
cl10_markers <- FindAllMarkers(cl10,
                               return.thresh = 0.01,
                               logfc.threshold = 0.5,
                               min.pct = 0.20,
                               verbose = F)

table(cl10_markers$cluster)

 0  1  2  3  4  5  6 
32 37 23 31 43 43 41 
cl10_markers %>%
  arrange(dplyr::desc(avg_log2FC)) %>%
  group_by(cluster) %>%
  dplyr::slice(1:5)
# A tibble: 35 × 7
# Groups:   cluster [7]
      p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  
      <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr> 
 1 1.09e-57       2.20 0.808 0.345  4.90e-55 0       MGP   
 2 6.47e-58       1.94 0.783 0.262  2.90e-55 0       DCN   
 3 3.78e-55       1.85 0.786 0.253  1.69e-52 0       CXCL14
 4 3.12e-18       1.60 0.363 0.128  1.40e-15 0       GPNMB 
 5 7.72e-22       1.55 0.491 0.205  3.46e-19 0       LUM   
 6 3.98e-66       3.70 0.541 0.064  1.78e-63 1       CCL14 
 7 1.33e-29       2.83 0.475 0.144  5.97e-27 1       POSTN 
 8 6.37e-66       2.53 0.847 0.258  2.85e-63 1       PLVAP 
 9 4.47e-11       1.96 0.224 0.07   2.00e- 8 1       ENTPD1
10 2.02e-29       1.92 0.557 0.171  9.05e-27 1       FLT1  
# ℹ 25 more rows
nos <- c()
fibro1 <- c() # FN1 high
fibro2 <- c(0, 4) # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c(3) # NOTCH3, RGS5
endo <- c(1, 2, 5) # PLVAP
myel <- c(6)

cl10$annotation <- ifelse(cl10$SCTsnn_res.0.8 %in% peri, "Peri_from_stromal",
                   ifelse(cl10$SCTsnn_res.0.8 %in% endo, "Endo_from_stromal",
                   ifelse(cl10$SCTsnn_res.0.8 %in% fibro2, "F2_from_stromal",
                   ifelse(cl10$SCTsnn_res.0.8 %in% myel, "Myel_from_stromal", NA))))

cl10$annotation <- paste0(cl10$SCTsnn_res.0.8, "_", cl10$annotation)

unique(cl10$annotation)
[1] "0_F2_from_stromal"   "5_Endo_from_stromal" "1_Endo_from_stromal"
[4] "4_F2_from_stromal"   "2_Endo_from_stromal" "6_Myel_from_stromal"
[7] "3_Peri_from_stromal"

Cluster 11

FeaturePlot(subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 11),
            features = c("FN1", "CEACAM6"),
            blend = T) &
  theme_bw()

reccl11 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 11)
reccl11 <- recluster_seurat(reccl11, umap_pcs = 10)

FeaturePlot(reccl11,
            features = c("DCN", "MARCO"),
            blend = T) &
  theme_bw()

FeaturePlot(reccl11,
            features = c("DCN", "MARCO")) &
  theme_bw()

ncol(reccl11)
[1] 1100
table(reccl11$SCTsnn_res.0.8)

  0   1   2   3   4   5   6   7   8 
240 205 199  97  87  85  75  69  43 
DimPlot(reccl11,
        group.by = "SCTsnn_res.0.8",
        reduction = "umap",
        label = T) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend() +
  ggtitle("Cluster 11")

reccl11$SCTsnn_res.0.8 <- factor(reccl11$SCTsnn_res.0.8,
                              levels = sort(unique(reccl11$SCTsnn_res.0.8)))
Idents(reccl11) <- reccl11$SCTsnn_res.0.8

DotPlot(reccl11,
        group.by = "SCTsnn_res.0.8",
        features = plot_features,
        cols = c("gray98", "tomato3")) +
  RotatedAxis()

nos <- c(3)
fibro1 <- c(0, 1, 2, 4, 7, 8) # FN1 high
fibro2 <- c() # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c() # NOTCH3, RGS5
endo <- c() # PLVAP
myel <- c(5, 6)

reccl11$annotation <- ifelse(reccl11$SCTsnn_res.0.8 %in% nos, "NOS_from_stromal",
                   ifelse(reccl11$SCTsnn_res.0.8 %in% fibro1, "F1_from_stromal",
                   ifelse(reccl11$SCTsnn_res.0.8 %in% myel, "Myel_from_stromal", NA)))

reccl11$annotation <- paste0(reccl11$SCTsnn_res.0.8, "_", reccl11$annotation)

unique(reccl11$annotation)
[1] "2_F1_from_stromal"   "0_F1_from_stromal"   "8_F1_from_stromal"  
[4] "3_NOS_from_stromal"  "1_F1_from_stromal"   "4_F1_from_stromal"  
[7] "7_F1_from_stromal"   "5_Myel_from_stromal" "6_Myel_from_stromal"

Cluster 13

FeaturePlot(subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 13),
            features = c("COL1A2", "APOE"),
            blend = T) &
  theme_bw()

reccl13 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 13)
reccl13 <- recluster_seurat(reccl13, umap_pcs = 10)

FeaturePlot(reccl13,
            features = c("COL1A2", "APOE"),
            blend = T) &
  theme_bw()

FeaturePlot(reccl13,
            features = c("COL1A2", "APOE")) &
  theme_bw()

ncol(reccl13)
[1] 910
table(reccl13$SCTsnn_res.0.8)

  0   1   2   3   4   5   6 
216 170 134 104 104  94  88 
DimPlot(reccl13,
        group.by = "SCTsnn_res.0.8",
        reduction = "umap",
        label = T) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend() +
  ggtitle("Cluster 13")

reccl13$SCTsnn_res.0.8 <- factor(reccl13$SCTsnn_res.0.8,
                              levels = sort(unique(reccl13$SCTsnn_res.0.8)))
Idents(reccl13) <- reccl13$SCTsnn_res.0.8

DotPlot(reccl13,
        group.by = "SCTsnn_res.0.8",
        features = plot_features,
        cols = c("gray98", "tomato3")) +
  RotatedAxis()

nos <- c()
fibro1 <- c(1, 2) # FN1 high
fibro2 <- c() # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c() # NOTCH3, RGS5
endo <- c() # PLVAP
myel <- c(0, 3, 4, 5, 6)

reccl13$annotation <- ifelse(reccl13$SCTsnn_res.0.8 %in% fibro1, "F1_from_stromal",
                   ifelse(reccl13$SCTsnn_res.0.8 %in% myel, "Myel_from_stromal", NA))

reccl13$annotation <- paste0(reccl13$SCTsnn_res.0.8, "_", reccl13$annotation)

unique(reccl13$annotation)
[1] "0_Myel_from_stromal" "2_F1_from_stromal"   "5_Myel_from_stromal"
[4] "3_Myel_from_stromal" "4_Myel_from_stromal" "1_F1_from_stromal"  
[7] "6_Myel_from_stromal"

Cluster 20

reccl20 <- subset(stromal_recluster_subset, subset = SCTsnn_res.0.8 == 20)
reccl20 <- recluster_seurat(reccl20, umap_pcs = 10)

FeaturePlot(reccl20,
            features = c("COL1A2", "PLVAP"),
            blend = T) &
  theme_bw()

FeaturePlot(reccl20,
            features = c("COL1A2", "PLVAP")) &
  theme_bw()

ncol(reccl20)
[1] 247
table(reccl20$SCTsnn_res.0.8)

 0  1  2  3 
80 78 48 41 
DimPlot(reccl20,
        group.by = "SCTsnn_res.0.8",
        reduction = "umap",
        label = T) +
  coord_fixed(ratio = 1) +
  theme_bw() +
  NoLegend() +
  ggtitle("Cluster 20")

reccl20$SCTsnn_res.0.8 <- factor(reccl20$SCTsnn_res.0.8,
                              levels = sort(unique(reccl20$SCTsnn_res.0.8)))
Idents(reccl20) <- reccl20$SCTsnn_res.0.8

DotPlot(reccl20,
        group.by = "SCTsnn_res.0.8",
        features = plot_features,
        cols = c("gray98", "tomato3")) +
  RotatedAxis()

nos <- c()
fibro1 <- c(0, 1, 3) # FN1 high
fibro2 <- c() # FN1 low, DCN, LUM (iCAF), COL1A2
peri <- c(2) # NOTCH3, RGS5
endo <- c() # PLVAP
myel <- c()

reccl20$annotation <- ifelse(reccl20$SCTsnn_res.0.8 %in% fibro1, "F1_from_stromal",
                   ifelse(reccl20$SCTsnn_res.0.8 %in% peri, "Peri_from_stromal", NA))

reccl20$annotation <- paste0(reccl20$SCTsnn_res.0.8, "_", reccl20$annotation)

Consolidating preliminary annotations

# Merging
seurat_list <- ls(pattern="reccl")
seurat_list <- str_sort(seurat_list, numeric = TRUE)
seurat_list <- do.call("list", mget(seurat_list))

#saveRDS(seurat_list, "/scratch/hnatri/PIPAC/seurat_list.Rds")
seurat_list <- readRDS("/scratch/hnatri/PIPAC/seurat_list.Rds")

seurat_list <- lapply(names(seurat_list), function(xx){
  #message(xx)
  cll <- xx
  xx <- seurat_list[[xx]]
  xx[["denoist_RNA"]] <- as(object = xx[["denoist_RNA"]], Class = "Assay5")
  DefaultAssay(xx) <- "denoist_RNA"
  xx[["RNA"]] <- NULL
  xx[["nucleus_RNA"]] <- NULL
  xx[["SCT"]] <- NULL
  xx$annotation <- paste0(cll, "_", xx$annotation)
  xx@meta.data <- xx@meta.data[,c("cell_id", "annotation")]
  
  xx
})

seurat_merged <- merge(x = seurat_list[[1]], y = seurat_list[2:3])

saveRDS(seurat_merged, "/scratch/hnatri/PIPAC/seurat_merged.rds")
#seurat_merged <- readRDS("/scratch/hnatri/PIPAC/seurat_merged.rds")

# Adding most granular annotations to the stromal reclustered subset object
stromal_recluster_subset$stromal_subcluster_annot <- mapvalues(x = colnames(stromal_recluster_subset),
                                                              from = colnames(seurat_merged),
                                                              to = seurat_merged$annotation)
stromal_recluster_subset@meta.data[-which(colnames(stromal_recluster_subset) %in% colnames(seurat_merged)),]$stromal_subcluster_annot <- NA

# Adding annotations to the complete stromal subset
stromal_ambig_annot <- read_sheet(metadata, sheet = "Stromal ambiguous recluster")
#stromal_ambig_annot <- stromal_ambig_annot[-which(is.na(stromal_ambig_annot$annotation)),]
#stromal_ambig_annot$annotation <- paste0(stromal_ambig_annot$SCTsnn_res.0.8, "_", stromal_ambig_annot$annotation)

stromal_seurat_data$stromal_subcluster_annot <- mapvalues(x = colnames(stromal_seurat_data),
                                                     from = colnames(stromal_recluster_subset),
                                                     to = stromal_recluster_subset$stromal_subcluster_annot)
stromal_seurat_data@meta.data[-which(colnames(stromal_seurat_data) %in% colnames(stromal_recluster_subset)),]$stromal_subcluster_annot <- NA

stromal_recluster_subset$stromal_ambig_annot <- mapvalues(x = stromal_recluster_subset$SCTsnn_res.0.8,
                                                         from = stromal_ambig_annot$SCTsnn_res.0.8,
                                                         to = as.character(stromal_ambig_annot$annotation))
#stromal_recluster_subset@meta.data[-which(stromal_recluster_subset$SCTsnn_res.0.8 %in% stromal_ambig_annot$SCTsnn_res.0.8),]$stromal_ambig_annot <- NA

# Adding stromal annotations to the full non-immune object
seurat_data$stromal_ambig_annot <- mapvalues(x = colnames(seurat_data),
                                             from = colnames(stromal_recluster_subset),
                                             to = as.character(stromal_recluster_subset$stromal_ambig_annot))
seurat_data@meta.data[-which(colnames(seurat_data) %in% colnames(stromal_recluster_subset)),]$stromal_ambig_annot <- NA

seurat_data$stromal_subcluster_annot <- mapvalues(x = colnames(seurat_data),
                                                  from = colnames(stromal_recluster_subset),
                                                  to = stromal_recluster_subset$stromal_subcluster_annot)
seurat_data@meta.data[-which(colnames(seurat_data) %in% colnames(seurat_merged)),]$stromal_subcluster_annot <- NA

# Adding labels for immune clusters in the whole stromal subset
stromal_seurat_data$stromal_annot <- ifelse(stromal_seurat_data$leiden_1.0 == 1, "Immune1_from_stromal",
                                     ifelse(stromal_seurat_data$leiden_1.0 == 11, "Immune2_from_stromal",
                                      ifelse(stromal_seurat_data$leiden_1.0 == 25, "Immune3_from_stromal", NA)))

seurat_data$stromal_annot <- mapvalues(x = colnames(seurat_data),
                                       from = colnames(stromal_seurat_data),
                                       to = stromal_seurat_data$stromal_annot)
seurat_data@meta.data[-which(colnames(seurat_data) %in% colnames(stromal_seurat_data)),]$stromal_annot <- NA

# Adding annotation for the whole non-immune subset
# first_pass_annot
seurat_data$nonimmune_annotation <- mapvalues(x = seurat_data$leiden_1.5,
                                       from = first_pass_annot$leiden_1.5,
                                       to = first_pass_annot$annotation)
seurat_data@meta.data[-which(seurat_data$leiden_1.5 %in% stromal_seurat_data$leiden_1.5),]$nonimmune_annotation <- NA

seurat_data$nonimmune_annotation <- as.character(seurat_data$nonimmune_annotation)

# Combining
#saveRDS(seurat_data, "/scratch/hnatri/PIPAC/nonimmune_annotated.rds")
#seurat_data <- readRDS("/scratch/hnatri/PIPAC/nonimmune_annotated.rds")

seurat_data$combined_annotation <- ifelse(is.na(seurat_data$stromal_subcluster_annot), seurat_data$stromal_ambig_annot, seurat_data$stromal_subcluster_annot)

seurat_data$combined_annotation <- ifelse(is.na(seurat_data$combined_annotation), as.character(seurat_data$stromal_annot), as.character(seurat_data$combined_annotation))

seurat_data$combined_annotation <- ifelse(is.na(seurat_data$combined_annotation), seurat_data$nonimmune_annotation, seurat_data$combined_annotation)

#seurat_data$combined_annotation <- ifelse(is.na(seurat_data$combined_annotation), seurat_data$annotation, #NA)

saveRDS(seurat_data, "/scratch/hnatri/PIPAC/nonimmune_annotated.rds")

# Rscript -e "rmarkdown::render('nonimmune_annotation_cell.Rmd')"
# mv *html /home/hnatri/PIPAC_spatial/docs/

sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ComplexHeatmap_2.18.0 viridis_0.6.3         viridisLite_0.4.2    
 [4] circlize_0.4.15       plyr_1.8.8            RColorBrewer_1.1-3   
 [7] googlesheets4_1.1.0   ggrepel_0.9.3         ggpubr_0.6.0         
[10] lubridate_1.9.2       forcats_1.0.0         stringr_1.5.0        
[13] purrr_1.0.1           readr_2.1.4           tidyr_1.3.0          
[16] tibble_3.2.1          ggplot2_3.4.2         tidyverse_2.0.0      
[19] SeuratDisk_0.0.0.9021 Seurat_5.0.1          SeuratObject_5.0.1   
[22] sp_1.6-1              dplyr_1.1.2           arrow_21.0.0.1       
[25] workflowr_1.7.1      

loaded via a namespace (and not attached):
  [1] fs_1.6.2                    matrixStats_1.0.0          
  [3] spatstat.sparse_3.0-1       bitops_1.0-7               
  [5] httr_1.4.6                  doParallel_1.0.17          
  [7] tools_4.3.0                 sctransform_0.4.1          
  [9] backports_1.4.1             utf8_1.2.3                 
 [11] R6_2.5.1                    lazyeval_0.2.2             
 [13] uwot_0.1.14                 GetoptLong_1.0.5           
 [15] withr_2.5.0                 gridExtra_2.3              
 [17] progressr_0.13.0            cli_3.6.1                  
 [19] Biobase_2.62.0              Cairo_1.6-0                
 [21] spatstat.explore_3.2-1      fastDummies_1.7.3          
 [23] labeling_0.4.2              sass_0.4.6                 
 [25] spatstat.data_3.0-1         ggridges_0.5.4             
 [27] pbapply_1.7-0               parallelly_1.36.0          
 [29] limma_3.58.1                rstudioapi_0.14            
 [31] generics_0.1.3              shape_1.4.6                
 [33] ica_1.0-3                   spatstat.random_3.1-5      
 [35] car_3.1-2                   Matrix_1.6-5               
 [37] ggbeeswarm_0.7.2            fansi_1.0.4                
 [39] S4Vectors_0.40.2            abind_1.4-5                
 [41] lifecycle_1.0.3             whisker_0.4.1              
 [43] yaml_2.3.7                  carData_3.0-5              
 [45] SummarizedExperiment_1.32.0 SparseArray_1.2.3          
 [47] Rtsne_0.16                  glmGamPoi_1.14.3           
 [49] promises_1.2.0.1            crayon_1.5.2               
 [51] miniUI_0.1.1.1              lattice_0.21-8             
 [53] cowplot_1.1.1               magick_2.7.4               
 [55] pillar_1.9.0                knitr_1.43                 
 [57] GenomicRanges_1.54.1        rjson_0.2.21               
 [59] future.apply_1.11.0         codetools_0.2-19           
 [61] leiden_0.4.3                glue_1.6.2                 
 [63] getPass_0.2-4               data.table_1.14.8          
 [65] vctrs_0.6.2                 png_0.1-8                  
 [67] spam_2.9-1                  cellranger_1.1.0           
 [69] gtable_0.3.3                assertthat_0.2.1           
 [71] cachem_1.0.8                xfun_0.39                  
 [73] S4Arrays_1.2.0              mime_0.12                  
 [75] survival_3.5-5              gargle_1.4.0               
 [77] iterators_1.0.14            statmod_1.5.0              
 [79] ellipsis_0.3.2              fitdistrplus_1.1-11        
 [81] ROCR_1.0-11                 nlme_3.1-162               
 [83] bit64_4.0.5                 RcppAnnoy_0.0.20           
 [85] GenomeInfoDb_1.38.5         rprojroot_2.0.3            
 [87] bslib_0.4.2                 irlba_2.3.5.1              
 [89] vipor_0.4.5                 KernSmooth_2.23-21         
 [91] colorspace_2.1-0            BiocGenerics_0.48.1        
 [93] ggrastr_1.0.2               tidyselect_1.2.0           
 [95] processx_3.8.1              bit_4.0.5                  
 [97] compiler_4.3.0              curl_5.0.0                 
 [99] git2r_0.32.0                hdf5r_1.3.8                
[101] DelayedArray_0.28.0         plotly_4.10.2              
[103] scales_1.2.1                lmtest_0.9-40              
[105] callr_3.7.3                 digest_0.6.31              
[107] goftest_1.2-3               spatstat.utils_3.0-3       
[109] presto_1.0.0                rmarkdown_2.22             
[111] XVector_0.42.0              htmltools_0.5.5            
[113] pkgconfig_2.0.3             sparseMatrixStats_1.14.0   
[115] MatrixGenerics_1.14.0       highr_0.10                 
[117] fastmap_1.1.1               rlang_1.1.1                
[119] GlobalOptions_0.1.2         htmlwidgets_1.6.2          
[121] DelayedMatrixStats_1.24.0   shiny_1.7.4                
[123] farver_2.1.1                jquerylib_0.1.4            
[125] zoo_1.8-12                  jsonlite_1.8.5             
[127] RCurl_1.98-1.12             magrittr_2.0.3             
[129] GenomeInfoDbData_1.2.11     dotCall64_1.0-2            
[131] patchwork_1.1.2             munsell_0.5.0              
[133] Rcpp_1.0.10                 reticulate_1.29            
[135] stringi_1.7.12              zlibbioc_1.48.0            
[137] MASS_7.3-60                 parallel_4.3.0             
[139] listenv_0.9.0               deldir_1.0-9               
[141] splines_4.3.0               tensor_1.5                 
[143] hms_1.1.3                   ps_1.7.5                   
[145] igraph_1.4.3                spatstat.geom_3.2-1        
[147] ggsignif_0.6.4              RcppHNSW_0.5.0             
[149] reshape2_1.4.4              stats4_4.3.0               
[151] evaluate_0.21               tzdb_0.4.0                 
[153] foreach_1.5.2               httpuv_1.6.11              
[155] RANN_2.6.1                  polyclip_1.10-4            
[157] future_1.32.0               clue_0.3-64                
[159] scattermore_1.2             broom_1.0.4                
[161] xtable_1.8-4                RSpectra_0.16-1            
[163] rstatix_0.7.2               later_1.3.1                
[165] googledrive_2.1.0           beeswarm_0.4.0             
[167] IRanges_2.36.0              cluster_2.1.4              
[169] timechange_0.2.0            globals_0.16.2